function v = deltamethod(fun,x,VCOV)

f = fun(x);

signx = sign(x);
signx(x==0) = 1;
h = sqrt(eps)*(signx.*max(abs(x),1));

G = zeros([numel(f) numel(x)]);
for i = 1:numel(x)
    x1 = x;
    x1(i) = x(i) + h(i);
    dx = x1(i) - x(i);
    f1 = fun(x1);
    G(:,i) = reshape((f1-f)./dx,[],1);
end

v = reshape(diag(G*VCOV*G'),size(f));
